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ABSTRACT: This analysis modifies the parsimonious specification of recently published total nitrogen (TN) and 
total phosphorus (TP) national-scale SPAtially Referenced Regressions On Watershed attributes models to allow 
each model coefficient to vary geographically among three major river basins of the conterminous United States. 
Regionalization of the national models reduces the standard errors in the prediction of TN and TP loads, 
expressed as a percentage of the predicted load, by about 6 and 7%. We develop and apply a method for combin- 
ing national-scale and regional-scale information to estimate a hybrid model that imposes cross-region con- 
straints that limit regional variation in model coefficients, effectively reducing the number of free model 
parameters as compared to a collection of independent regional models. The hybrid TN and TP regional models 
have improved model fit relative to the respective national models, reducing the standard error in the prediction 
of loads, expressed as a percentage of load, by about 5 and 4%. Only 19% of the TN hybrid model coefficients 
and just 2% of the TP hybrid model coefficients show evidence of substantial regional specificity (more than 
±100% deviation from the national model estimate). The hybrid models have much greater precision in the esti- 
mated coefficients than do the unconstrained regional models, demonstrating the efficacy of pooling information 
across regions to improve regional models. 
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INTRODUCTION 



River Basin and Gulf of Mexico) where the develop- 
ment of effective solutions is complicated by the 
diversity of nutrient contributions from upstream 
sources and catchments (New York State Department 
of Environmental Conservation and Connecticut 



Water managers increasingly have a need for 
model-based information to address eutrophication 
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the watershed response at these large scales places 
multiple demands on the water-quality model. The 
model must provide an accurate assessment of water- 
quality conditions for locations that are deficient in 
monitoring data; the model must make an accurate 
determination of the relative importance of an exten- 
sive list of nutrient sources contributing to a particu- 
lar water-quality impairment, many of these sources 
being physically distant from the impaired area; and 
the model must adequately describe how the aquatic 
system will respond to hypothetical changes in nutri- 
ent sources or factors affecting their transport. 

The U.S. Geological Survey's (USGS) SPAtially 
Referenced Regressions On Watershed attributes 
(SPARROW) model - a hybrid statistical and process- 
based mass balance model - was developed to explain 
spatial variability in mean water-quality conditions 
and predict source contributions in the streams of 
large, heterogeneous basins (Schwarz etal., 2006; 
Preston et al., 2009) (see also the brief description of 
the SPARROW methodology contained in the Sup- 
porting Information for this article). SPARROW has 
been used previously to assess nutrient loadings in 
large watersheds in the United States (U.S.) (e.g., 
Smith etal., 1997; Alexander etal., 2008) and other 
countries (Alexander et al., 2002). An important ques- 
tion in applying SPARROW and other water-quality 
models in large watersheds is how the governing 
material transport equations should differ spatially 
in their functional forms and/or coefficients to 
account for the effects of landscape heterogeneity and 
spatial scaling on nutrient supply and transport pro- 
cesses. Despite advances in watershed science to 
address this question, uncertainties remain over the 
nature of these effects and how they should be repre- 
sented in watershed models (e.g., McDonnell etal., 
2007; Kirchner, 2009). Existing national-scale 
SPARROW nutrient models are estimated with the 
assumption that the coefficients of the model that 
describe the marginal effect of water-quality explana- 
tory variables on stream nutrient loads (e.g., the 
loading response from an incremental change in soil 
permeability; the change in the fraction of nutrient 
removed in streams from a unit change in water tra- 
vel time) are the same for all locations nationally 
(Alexander et al., 2008). This assumption insures that 
the widest ranges of environmental conditions are 
used to estimate the effects of individual model pro- 
cesses, and imparts a parsimonious form that facili- 
tates the interpretation of the estimated coefficients. 
The assumption may be overly simplistic, however, if 
there exist latent processes, related to landscape het- 
erogeneity or process interactions, for which surro- 
gate measures of their effects are not available in 
geospatial datasets. Under these conditions, large- 
scale model development could lead to biased predic- 



tions in particular regions and under certain manage- 
ment or forecasting scenarios. 

In recognition of potential regional variation in 
water-quality model processes, the USGS National 
Water Quality Assessment (NAWQA) Program has 
developed a set of SPARROW regional major river 
basin (MRB) nutrient models, described in this Fea- 
tured Collection (for a summary of these models, see 
Preston et al., this issue). The estimation of these 
regional models has benefited from a large assem- 
blage of water-quality data, representing an approxi- 
mate order of magnitude increase in the number of 
monitoring sites, across all regions, compared to the 
number of sites from which data were available to 
construct existing national-scale SPARROW nutrient 
models (Smith etal., 1997; Alexander etal., 2008). 
The inherent flexibility afforded by independent 
regional models allows the predictions of current 
water-quality conditions derived from these models to 
have less bias and greater precision than those of 
national models; thus, from the standpoint of water- 
quality assessment, the regional models out-perform 
the national models. 

Regionalization of models, however, also poses a 
potential negative consequence: an individual regional 
model likely encompasses a narrower range in basin 
attributes than does a national model. For example, a 
climate variable typically exhibits less variation within 
a region than across the entire nation. In a model, an 
attribute's effect on water quality is mediated by a 
coefficient, the estimation of which is made more pre- 
cise the greater the variation of the attribute in the 
sample of data used to estimate the model. The 
reduced variation of an attribute within a confined 
region thereby reduces the precision of the attendant 
coefficient. Thus, users of a regional model (e.g., water 
resource managers) face a fundamental tradeoff in 
model accuracy between the bias and precision of the 
model predictions. On the one hand, if a model coeffi- 
cient's true value varies across regions, the loss in pre- 
cision due to less within-region variability in the 
associated basin attribute may be an acceptable out- 
come for obtaining an unbiased coefficient estimate. 
On the other hand, if the true coefficient does not vary 
across regions, the reduced precision implies a loss in 
accuracy of the estimated coefficient, a loss that possi- 
bly could be ameliorated by expanding the spatial 
extent of the model to include a wider range of basin 
attributes. It remains an open question, therefore, 
whether the advantage of added flexibility afforded by 
individual regional models, which facilitates a less 
biased assessment of existing water-quality conditions, 
outweighs the loss of precision, which potentially 
diminishes a water manager's ability to use the model 
for predicting changes in water quality brought about 
by changes in its causative factors. 
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This article has two major objectives: the first 
objective is to evaluate regional variation in the 
nutrient model coefficients using a "fixed-effects" 
approach (Judge et al., 1985) to determine whether a 
single national model should be rejected in favor of a 
more regional-based approach. The "fixed-effects" 
approach to regionalization requires specifying multi- 
ple sets of coefficients in the national model, a sepa- 
rate set for each region of the analysis, with each 
regional set consisting of all the process coefficients 
included in the national model. This approach is a 
generalization of the fixed effects approach commonly 
used in longitudinal studies, which typically limit the 
fixed effects specification to an intercept term. Func- 
tionally, with fixed-effects specified for all process 
coefficients, the approach is equivalent to estimating 
independent regional models, the manner of regional- 
ization implemented by the SPARROW MRB models. 
The fixed-effects method has been used, to varying 
degrees, in numerous previous SPARROW modeling 
studies (for examples of partial approaches, whereby 
fixed effects are specified only for a subset of the pro- 
cess coefficients, see Smith etal., 2003; Hoos and 
McMahon, 2009). An advantage of the approach is its 
flexibility — it places no restrictions on the pattern of 
variation taken by the regional process coefficients. 
The chief disadvantage, as mentioned above, is that 
the precision of the coefficient estimates of the inde- 
pendently estimated models may be compromised if 
there is too little within-region variation in basin 
attributes. 

The second major objective of this article is to eval- 
uate the potential for improving the precision of the 
regional coefficient estimates. As suggested above, if 
it is known that regional coefficients are similar, it 
may be possible to improve precision by pooling infor- 
mation across regions. In the context of the fixed- 
effects approach, one way to do this is to estimate the 
model with cross-region constraints placed on the 
coefficients. The resulting model, which we call a 
hybrid model, effectively pools all the region-specific 
information, allowing the estimation of process coeffi- 
cients that are not statistically different across 
regions to benefit from the greater variation in basin 
attributes afforded by a national-scale analysis. The 
constrained estimates of the fixed-effects coefficients 
retain some degree of regional variation, with the 
extent of variation depending on the particular cross- 
region constraints imposed as part of the specification 
of the hybrid model. 

A distinct advantage of the hybrid model, from the 
standpoint of the analysis described here, is that it 
encompasses both the regional fixed-effects model 
and the national model. At one extreme, if there 
are no constraints, the hybrid model is identical to 
the regional fixed-effects model, the model that is 



functionally equivalent to specifying independent 
regional models. At the other extreme, if all possible 
independent constraints are applied, the hybrid 
model becomes identical to the national model. Thus, 
by partially constraining the regional fixed-effects 
model, the hybrid model provides a framework for 
evaluating the tradeoff between achieving prediction 
accuracy and improving coefficient precision in a 
regional context. 

The national models serving as the subject of the 
fixed-effects analysis are the previously published 
SPARROW models for total nitrogen (TN) and total 
phosphorus (TP), estimated for the base year 1992 
(Alexander etal., 2008). As our interest is to isolate 
the effects of regionalization, yet retain as much as 
possible the specifications of the existing national 
models, this analysis uses data for the same set of 
monitoring stations that were used by Alexander 
et al. (2008). This collection of stations is too small, 
however, to obtain reasonable regional estimates for 
the number of regions used in the NAWQA MRB 
regional studies, which are based on many more sta- 
tions. To compensate, this study adopts a rather 
coarse partitioning of the nation into only three 
regions: the East, Northwest, and Southwest regions 
of the conterminous U.S. (see Figure 1). This coarse 
regionalization makes it unlikely that the full range 
of regional effects present in the U.S. are properly 
accounted for, and supporting evidence of this short- 
coming is presented. The limited size of the dataset 
used in this analysis also makes it likely that some of 
the regional variation exhibited by the hybrid model 
coefficient estimates is spurious. For these reasons, it 
is prudent to view the hybrid model coefficient esti- 
mates presented here as preliminary. The subsequent 
production of "management-ready" coefficient esti- 
mates should be possible by applying the hybrid 
method to the expansive set of data incorporated in 
the NAWQA MRB studies, the analysis of which will 
likely also benefit from altering the specification of 
the models beyond the existing national SPARROW 
models. Despite the shortcomings, we believe the 
data are sufficient to demonstrate the potential util- 
ity of moving beyond modeling frameworks that are 
either exclusively national or strictly regionally 
independent. 

The following analysis begins with a brief descrip- 
tion of the recently published national nutrient 
SPARROW models (Alexander et al, 2008), developed 
for the l:500,000-scale Reach File 1 network (Nolan 
et al., 2002). An initial evaluation of the national 
nutrient models is undertaken to determine if regio- 
nal fixed effects are present in three major regions of 
the conterminous U.S., with attention given to 
whether the regional effects can be associated with 
specific groups of transport processes that are an 
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A. Total Nitrogen (TN) 




▼ 1.0 to 2.5 



B.Total Phosphorus (TP) 




T 1.5 to 4.3 



FIGURE 1. The Residuals Estimated From the 1992 National 
SPARROW Models for (A) Total Nitrogen and (B) Total Phosphorus 
(Alexander et al., 2008). The residuals for total nitrogen are 
weighted to account for uncertainty in the monitored load esti- 
mates (Alexander et al., 2008, p. S-12). A negative residual implies 
overprediction of the SPARROW model. The number of monitoring 
sites in the region is given in parentheses. 

intrinsic feature of a SPARROW model specification. 
A subsequent section describes a practical approach 
for identifying constraints to be applied to the regio- 
nal fixed-effects model to estimate the hybrid model. 
A section describing the results of the estimation of 
the hybrid model follows, with emphasis placed on 
demonstrating the effects of agglomerating informa- 
tion across regions on the precision of coefficient esti- 
mates and model predictions. The concluding sections 
contain a discussion of regionalization methodology 
in the context of other approaches found in the litera- 
ture and a summary of the findings of this study. All 
mathematical derivations used in the analysis are 
described in an Appendix and in the Supporting 
Information for this article. Additionally, the Sup- 
porting Information includes a brief description of the 
SPARROW model and a detailed presentation of 
results from the regional fixed-effects and hybrid 
models. 



THE NATIONAL SPARROW NUTRIENT MODELS 



The specific TN and TP SPARROW models serving 
as the subject of the regionalization analysis are 
described in Alexander et al. (2008). Both models 
adopt 1992 as the base year for the analysis, and use 
the spatial infrastructure provided by the recent ver- 
sion of the enhanced Reach File 1 (E2RF1) stream 
network (Nolan et al., 2002), a l:500,000-scale net- 
work consisting of approximately 65,000 reach seg- 
ments in the conterminous U.S., including digital 
elevation-delineated catchments for each reach. 

The SPARROW modeling approach groups the 
determinants of water quality into three types: source 
variables that identify specific sources of contami- 
nant, land-to-water variables that mediate the 
amount of contaminant source that is transported to 
streams included in the digital stream network, and 
aquatic removal rate variables that determine the 
amount of contaminant that is transported through 
streams in the digital network (a detailed description 
of the SPARROW methodology is contained in the 
Supporting Information; see Schwarz et al., 2006). 
The nutrient sources used in the national nutrient 
models include long-term mean atmospheric deposi- 
tion, detrended to 1992 (TN model only); 1990 popula- 
tion (an urban source surrogate); 1992 National Land 
Cover Data estimates of forested land, barren land, 
and shrub land; and 1992 manure and fertilizer nitro- 
gen and phosphorus loading to agricultural lands sep- 
arately identified as growing corn and soybeans, 
alfalfa, wheat (TN model only), and other crops. The 
land-to-water delivery variables include soil perme- 
ability, drainage density (TN model only), mean tem- 
perature (TN model only), mean precipitation, 
specific catchment area (the catchment average of the 
area upslope of each grid cell within the catchment, 
divided by the average cell width; see the Supporting 
Information of Alexander et al., 2008), percent of area 
artificially drained, and mean catchment slope (TP 
model only). The specified models include proportion- 
ate loss processes in reach network stream and reser- 
voir segments. Both the stream and reservoir 
processes are based on empirically estimated mean 
settling velocities, a conceptual mechanism affecting 
the mean rate at which contaminants are removed 
from the aquatic environment (Schwarz et al., 2006). 

The models are estimated using long-term mean 
annual loads computed from monitoring data com- 
piled at 425 stations located on the E2RF1 stream 
network. The mean load estimates are based on 
water-quality observations collected over the period 
1975-1995, and daily streamflow data from the period 
1975-2000. The mean load estimates are detrended 
to the base year 1992 using methods described in 
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Schwarz et al. (2006). Each observation in the TN 
model is weighted to account for variations in the 
precision of the mean load estimate (procedure 
described in Alexander et al., 2008, p. S-12); all obser- 
vations in the TP model receive an equal weight. The 
land-to-water delivery variables are expressed as dif- 
ferences from their mean, a transformation that aids 
in the interpretation of the source coefficients, and 
will enable the comparison of the national models 
with the regional models described in the next sec- 
tion. Based on physical considerations, all source and 
aquatic delivery coefficients are constrained to be 
nonnegative. A complete description of the TN and 
TP models, including references to data sources and 
details on the methods used to create the dependent 
and explanatory variables, is provided in Alexander 
et al. (2008). 



EVALUATION OF NATIONAL NUTRIENT 
MODELS FOR REGIONAL EFFECTS 



A coarse assessment of regional effects in the 
national nutrient models can be made by visual 
inspection of the national model residuals (Alexander 
et al., 2008), displayed in Figure 1. The residuals are 
computed in natural logarithm space, the space in 
which the SPARROW models are estimated, and for 
TN are weighted to account for variations in the pre- 
cision of the load estimate (the observational weights 
noted above and described in Alexander et al. , 2008). 
A negative residual implies the SPARROW model 
overpredicts measured load. The most notable fea- 
ture, evident with both the TN and TP residuals, is 
the near absence of large residuals in the East. The 
SPARROW models are considerably less precise in 
arid environments, notably along the large climate 
gradient identified with the 100th meridian where 
large positive and negative errors for both nutrients 
are prominent. Regional bias is also evident, particu- 
larly for TP, which shows a preponderance of large 
positive residuals in the Northwest region and large 
negative residuals in the Southwest region. Although 
much less noticeable, the pattern of regional bias is 
similar for TN, with a prevalence of positive residuals 
in the Northwest region and negative residuals in the 
Southwest region. Also discernable are residual pat- 
terns within the regions. In the East region, TN 
shows a clustering of negative residuals along the 
South Atlantic coast (a pattern that was previously 
noted by Hoos and McMahon, 2009), and a grouping 
of positive residuals in Illinois/Iowa. A prevalence of 
positive TN residuals is also observed along the Paci- 
fic coast, in both the Northwest and Southwest 



regions. Wi thin-region patterns are more apparent 
with TP. A cluster of negative residuals is evident in 
the Lower Great Lakes (within the East region), 
along the Gulf coast (in both the East and Southwest 
regions), and the area extending across central Texas, 
eastern New Mexico, and western Oklahoma (South- 
west region). Clusters of positive TP residuals are 
observed in the East region in southern Florida, 
through much of the Midwest (in both the East and 
Northwest regions), and the eastern portions of the 
Southwest region. The existence of these intra-regio- 
nal residual patterns underscores the coarseness of 
the three-region characterization of regional specific- 
ity adopted in this study. 

A statistical assessment of regional patterns can be 
made by evaluating the national model residuals for 
spatial correlation. One such statistic that has been 
used for this purpose is Moran's / (Cliff and Ord, 
1973), evaluated here using spatial influence factors 
applied to the residuals, the factors being the inverse 
of the distance between paired sites, a weighting 
scheme that is more sensitive to large-scale 
regional correlation than to local correlation (Hoos 
and McMahon, 2009). The Moran's J has a statistical 
significance of <0.001 for both the TN and TP models, 
providing strong evidence for regional-scale spatial 
correlation. 

To quantify the regional patterns exhibited by the 
residuals in Figure 1, we developed a fixed-effects 
regionalization of the national SPARROW models by 
re-specifying the models to include three region-spe- 
cific coefficients for each process coefficient. As 
explained in the Introduction, this approach to region- 
alization is functionally equivalent to estimating an 
independent SPARROW model for each region (as long 
as the regional basins are independent, or if not, as is 
the case with the East region, part of which is down- 
stream from the Missouri basin component of the 
Northwest region, as long as region boundaries are 
located at a monitoring station). A formal statistical 
evaluation of the effects of regionalization on predic- 
tion accuracy was undertaken by comparing measures 
of fit between the national and regionalized nutrient 
models. Table 1 presents the difference between the 
root mean squared error (RMSE) of the national nutri- 
ent model and the alternative models with varying 
degrees of regionalization specified in the model pro- 
cess coefficients. The statistical significance of these 
differences is evaluated using standard .F-tests, with 
the national models representing the null hypothesis 
of no regionalization. Process coefficients are grouped 
according to landscape and aquatic properties, with 
landscape coefficients corresponding to the com- 
bined source and land-to-water delivery parameters, 
and aquatic process coefficients corresponding to 
the stream and reservoir removal-rate parameters. 
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TABLE 1. The Reduction in the Root Mean Squared Error (RMSE) of Regional Models Incorporating Varying 
Degrees of Region-Specific Fixed Effects for Land and Aquatic Processes, as Compared to the 1992 National 
SPARROW Models for Total Nitrogen (TN) and Total Phosphorus (TP) (Alexander el al, 2008). 



Alternative Model Specification TN TP 

Reduction in 



Regional 
Effects 


Process Coefficients 
for Which Region-Specific 
Fixed Effects Are 
Incorporated 


Reduction in RMSE 
Compared to National 
Model Without Regional 
Fixed Effects 1 


Number of 
Fixed-Effects 
Coefficients* 


RMSE Compared 
to National Model 
Without Regional 
Fixed Effects 1 


Number of 
Fixed-Effects 
Coefficients* 


All regions 


Land and aquatic 


0.055** 


54 


0.073** 


45 




Land only 


0.046** 


50 


0.058** 


41 




Aquatic only 


0.008* 


22 


0.033** 


19 


East region only 


Land and aquatic 


0.025** 


36 


0.039** 


30 




Land only 


0.026** 


34 


0.037** 


28 




Aquatic only 


<0.001 


20 


<0.001 


17 


Northwest 


Land and aquatic 


0.026** 


36 


0.044** 


30 


region only 


Land only 


0.018** 


34 


0.024** 


28 




Aquatic only 


0.007* 


20 


0.034** 


17 


Southwest 


Land and aquatic 


0.035** 


36 


0.053** 


30 


region only 


Land only 


0.032** 


34 


0.041** 


28 




Aquatic only 


0.007* 


20 


0.029** 


17 



Notes: Statistical significance of a reduction is determined using an F-test, with **denoting significance at the 0.001 level, *denoting 
significance at the 0.01 level; RMSE reductions without asterisks are not significant at the 0.05 level. 

Calculated as the national model RMSE minus the regional fixed-effects model RMSE. Reduction in RMSE of the natural log residuals can 
be interpreted as the approximate reduction in the standard error of a prediction of load for a reach, expressed as a share of the predicted 
load in the reach. 

'"The national model for TN has 18 coefficients. 
§ The national model for TP has 15 coefficients. 



Separate regional effects are specified for the identified 
process coefficient groups for either "all regions," 
meaning that each region has a unique effect for the 
stated process coefficient group, or for an individual 
region, meaning that only the identified region has a 
unique effect for the stated group; all other regions 
and unstated process coefficient groups are restricted 
so as to have a common effect. 

The results of these tests, as presented in Table 1, 
show that the national model specification, for either 
TN or TP, is strongly rejected in favor of fully region- 
alized models involving both land and water process 
coefficients jointly, or with land and water processes 
coefficients regionalized separately. As the reduction 
in the RMSE, multiplied by 100, can be interpreted 
as the approximate reduction in the standard error of 
a prediction of load for any given reach, expressed as 
a percent of predicted load, the results imply a fully 
regionalized model improves predictions by 6% for 
TN and by 7% for TP. Comparisons with models hav- 
ing only single region effects for combined land and 
water processes (the first line of each group in 
Table 1, excluding the first group) show rejections of 
the national model at a 0.001 significance level for all 
cases, with reductions in RMSE ranging from 3 to 4% 
for TN and 4 to 5% for TP. 

The comparative reductions in RMSE between TN 
and TP are in general agreement with Figure 1, 



which shows a more pronounced regional pattern for 
TP than for TN. The meager improvements in RMSE 
obtained with the regional fixed-effects model, for 
both TN and TP, is partly a consequence of the 
coarse, three-region delineation of the national model 
used in this analysis. As shown in Figure 1, the 
regions lack uniformity with respect to the residuals 
and it may be possible to explain more variation by 
modifying the delineation of regions to isolate groups 
of residuals having similar values. However, achiev- 
ing a large reduction in RMSE by partitioning the 
U.S. into more regions, delineated on the basis of the 
residuals depicted in Figure 1, is not guaranteed 
because SPARROW does not include a constant inter- 
cept (Schwarz et al. , 2006); any adjustment in a coef- 
ficient to correct regional bias must interact with a 
spatially varying land or aquatic causal variable, the 
spatial homogeneity of which is not necessarily coin- 
cident with that of the residuals. 

The significance of comparisons between the 
national models and individual region fixed-effects 
models with transport process groups considered sep- 
arately is also shown in Table 1 (the second and third 
lines of each group). These results show that land 
transport processes generally exhibit more region- 
specific variation than aquatic processes. The reduc- 
tion in error with regionalization applied to land 
processes alone across all regions is greater than the 
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reduction for aquatic process regionalization, with 
this difference being relatively smaller for TP (for 
TN, a 5% reduction due to land processes vs. a 1% 
reduction due to aquatic processes; for TP, a 6% 
reduction for land processes vs. a 3% reduction for 
aquatic processes). Land processes are statistically 
significant, at the 0.05 level, in all individual regions 
for both TN and TP. Conversely, TN and TP aquatic 
processes are statistically significant at the 0.05 level 
only for the individual Northwest and Southwest 
regions. Nutrient land processes for individual 
regions generally have a larger effect on RMSE than 
aquatic processes in all regions, except for TP in the 
Northwest region. 

The statistical significance of the Moran's J statis- 
tic for the residuals based on the complete regionali- 
zation of model processes is 0.144 for TN and 0.023 
for TP (see Table 2), indicating that, at least for TN, 
the three-region delineation of the national model 
has successfully removed spatial correlation that is 
sensitive to an inverse distance spatial influence fac- 
tor. The more pronounced intra-regional patterns in 
the TP residuals, noted in Figure 1, suggest that 
removal of spatial correlation requires a partitioning 
of the U.S. into more than three regions, although 
such an approach may not be successful for the 



reason explained above that each SPARROW coeffi- 
cient must interact with a spatially varying land or 
aquatic causal variable. 



DESCRIPTION OF THE HYBRID MODEL 



The principal conclusion to be drawn from the pre- 
vious section is that regional variation in the process 
coefficients is empirically confirmed: the model fits of 
the regional fixed-effects nutrient models are signifi- 
cantly improved, statistically, relative to those of the 
national models, although the quantitative improve- 
ments are not substantial. The emphasis now turns 
to regionalizing the national model in such a way 
that regional coefficients are estimated as precisely 
as possible without a large compromise in model pre- 
cision. This section describes a method for regionaliz- 
ing a national model using the fixed-effects approach 
by including in the specification a set of cross-region 
constraints on the fixed-effects coefficients. The 
resulting model, called a hybrid model, is intended to 
have more precise coefficient estimates than the 
regional fixed-effects model, but also reflect regional 



TABLE 2. Summary Statistics for the National, Regional Fixed-Effects, 
and Hybrid SPARROW Models for Total Nitrogen (TN) and Total Phosphorus (TP). 







TN 






TP 








Regional 






Regional 






National 


Fixed-Effects 


Hybrid 


National 


Fixed-Effects 


Hybrid 




Model 


Model 


Model 


Model 


Model 


Model 


Number of observations 


425 


425 


425 


425 


425 


425 


Number of coefficients 


18 


54 


54 


15 


45 


45 


Number of coefficients with 


0 


5 


4 


0 


1 


1 


binding physical constraints 














Number of cross-region constraints 


0 


0 


22 


0 


0 


28 


Error degrees of freedom 


407 


376 


397 


410 


381 


409 


Percent of coefficients with p < 0.05* 


72 


57 


76 


80 


41 


80 


Average coefficient of variation ' 


0.403 


0.891 


0.549 


0.369 


2.416 


0.389 


Sum of squared errors (SSE) 


124.3 


93.3 


101.3 


238.1 


180.9 


215.8 


Mean squared error (MSE) 


0.305 


0.248 


0.255 


0.581 


0.475 


0.528 


Root mean squared error (RMSE) 


0.553 


0.498 


0.505 


0.762 


0.689 


0.726 


R 2 


0.933 


0.953 


0.945 


0.870 


0.902 


0.883 


Adjusted R 2 


0.930 


0.947 


0.942 


0.866 


0.891 


0.878 


Yield R 2% 


0.866 


0.906 


0.891 


0.684 


0.760 


0.714 


Median absolute pet. 


42.3 


40.2 


38.4 


67.1 


52.6 


59.1 


prediction error 8 














Significance of Moran's I 


<0.001 


0.144 


0.026 


<0.001 


0.023 


<0.001 



*For the calculation of the percent of coefficients with p < 0.05, a one-tailed p-value is used for coefficients subject to physical bounds (all 
source and aquatic removal coefficients); otherwise, a two-tailed p-value is used. 

'The standard error of a coefficient estimate divided by its estimated value, averaged over all coefficients in the model. 

"'Yield R 2 adjusts the R 2 to account for the area of each observation's upstream basin, thereby removing from R 2 the inflating effects of a wide 
range of basin scales (Schwarz et al., 2006). 

^Median absolute error in prediction, expressed as a percent of the monitored load, for 678 TN and 865 TP stations not used in the 
estimation of the national, regional fixed-effects, or hybrid models. 
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variation in targeted coefficients sufficient to exhibit 
improved model fit as compared to the national 
model. 

Figure 2 summarizes the relationships linking the 
regional fixed-effects model with the hybrid model. 
There are three principal steps required to derive the 
hybrid model: the identification of a set of potential 
cross-region linear constraints, a statistical evalua- 
tion to determine which of these constraints are 
valid, and the estimation of the fixed-effects model 
with the valid constraints imposed. As shown in Fig- 
ure 2, the covariances of the regional fixed-effects 
model coefficient estimates are used to define a set of 
weights, these weights forming the constant terms of 
the potential linear constraints. The weights are also 
used, in conjunction with the regional model fixed- 
effects coefficient estimates, to form a set of comple- 
ment-region weighted-average estimates, a separate 
average for each fixed-effect coefficient. A stepwise 
procedure is used to compare each fixed-effect coeffi- 
cient with its complement-region weighted average 
counterpart, and the outcome of that comparison 
determines which of the potential linear constraints 
are applied in the estimation of the hybrid model. 
The following discussion explains in detail the moti- 
vation underlying each of these relationships. 

The imposition of constraints on the regional fixed- 
effects coefficients begs a question: What specific form 



Regional fixed-effects model 



Region-specific coefficients 



Region-specific covariances 



Complement-region weighted- 
average coefficient estimates 



Region-specific 
weights 



Stepwise procedure used to 

evaluate the statistical 
significance of differentials 
representing the differences 

between step-updated 
individual regional models 
and the complement-region 
weighted average estimates 



Region-specific weights 
define the constants of a set 

of linear constraints the 
empirical validity of which are 
evaluated using the stepwise 
procedure 



Hybrid model 



FIGURE 2. Flow Diagram Describing the Intermediate 
Steps Applied to the Results of the Regional Fixed-Effects 
Model to Obtain the Hybrid Model. 



should the set of constraints take? Many alternatives 
are possible, perhaps the most general of which is the 
set of all possible "pair-equivalence" constraints that 
equilibrate the process coefficient in one region to a 
like-process coefficient in another region (e.g., the 
atmospheric deposition coefficient for the East region 
equals the atmospheric deposition coefficient for the 
Northwest region, etc.). If there are R regions, and if 
the national model has K process coefficients, this 
general set consists of N c = K x R x (R - l)/2 differ- 
ent constraints (the K coefficients of region 1 equated 
with the K coefficients of each of the R - 1 other 
regions, plus the K coefficients of region 2 equated 
with the K coefficients of each of the R - 2 other 
regions excluding region 1, etc.), although not all of 
these constraints are independent. If all constraints 
are binding, meaning that all process coefficients in 
any region are equal to the corresponding like-process 
coefficients in any other region, the resulting model 
is equivalent to the national model. The practical dif- 
ficulty with this general constraint set concerns the 
computation time required to evaluate all of the con- 
straints. There are unique combinations of the iV c 
possible constraints, and the evaluation of each com- 
bination requires a separate estimation of the Kx R 
coefficients of a constrained nonlinear SPARROW 
fixed-effects model. 

An alternative approach is to adopt the perspective 
of a regional water-quality manager tasked with the 
problem of specifying a regional SPARROW model. 
We suppose that the manager is familiar with the 
regional fixed-effects SPARROW model estimated for 
the entire nation and has an interest in using the 
results from this model to improve the precision of 
some of the coefficients in the subject-region's model. 
The only restriction we place on the manager is that 
only the model for the subject region can be re- 
estimated; the regional fixed-effects coefficient esti- 
mates for the complement regions cannot be altered 
(ultimately, this restriction is removed in the estima- 
tion of the hybrid model, which involves the re-esti- 
mation of all fixed-effect coefficients across all 
regions). The proposed approach is to compare the 
estimate of a process coefficient in the subject region 
to a weighted average of the fixed-effect coefficients 
for the complement regions. If the subject-region's 
estimate is not statistically different from the 
weighted-average estimate, the weighted-average esti- 
mate is substituted for the subject-region's estimate, 
effectively constraining the subject-region's coefficient 
to equal the weighted-average estimate; otherwise, 
the subject-region's coefficient remains a free parame- 
ter to be estimated in the subject-region model. 

Of practical concern is the determination of the 
weights used to form the weighted average of the 
complement-region fixed-effect coefficient estimates. 
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In the Appendix it is shown that there exists a 
unique set of weights, derived from the covariance 
matrices of the complement-region fixed-effect coeffi- 
cient estimates, that minimizes the variance of the 
difference between the subject-region coefficient 
estimate and the associated complement-region 
weighted-average estimate. These weights have a 
number of interesting features and statistical impli- 
cations. First, the minimization of the variance of the 
differential implies that the power of the statistical 
test comparing the subject-region coefficient estimate 
to the corresponding weighted-average estimate is 
maximized, thereby improving the probability of 
detecting a statistically significant difference in the 
subject-region's coefficient estimate. Thus, the speci- 
fied weights serve to minimize the number of con- 
straints imposed in both the subject-region model 
and, as explained below, the hybrid model. 

Second, because the covariance matrices typically 
include nonzero covariances for all pairings of the 
process coefficients, the optimal weights are nonzero 
for all complement-region coefficients, not just for 
like-process coefficients in these regions. Moreover, 
although the weights sum to one and are necessarily 
positive for like-process coefficients in the comple- 
ment regions, the weights may be either positive or 
negative for dissimilar-process coefficients, depending 
in a complicated way on the signs and values of the 
individual elements of the covariance matrices. Thus, 
for example, in the determination of the complement- 
region weighted average for the atmospheric deposi- 
tion coefficient in the East region, the weights will be 
positive for the atmospheric deposition coefficients in 
the Northwest and Southwest regions, but could be 
either positive or negative for, say, the reservoir 
removal rate coefficients in those regions. As shown 
in the Appendix, the weights for dissimilar-process 
coefficients are all zero only if the covariance matri- 
ces across all complement regions are identical up to 
a scaling factor equal to the number of water-quality 
monitoring sites in the region. Under such conditions, 
the optimal weights for like-process coefficients are 
equal to the share of monitoring sites in each region. 

A third interesting feature of the optimal weights, 
also derived in the Appendix, concerns the intuitive 
interpretation of the weighted-average estimates if 
the underlying regional model is linear in the coeffi- 
cients. In this case, the weighted-average estimates 
reduce identically to the estimates that would result 
if all like-process coefficients in the complement 
regions are constrained to have the same value, these 
values being estimated by weighted least squares 
with weights equal to the inverse of the variance of 
the model residuals in each region. Thus, in a linear 
model, the complement-region weighted-average esti- 
mate essentially loses all regional specificity, and the 



comparison of the subject-region coefficient to the 
complement-regions' weighted average is simply a 
comparison of two regional estimates - the subject- 
region estimate and the corresponding estimate for 
the agglomerated complement regions. This equiva- 
lence fails to hold for a nonlinear model such as 
SPARROW because the covariance matrices are 
themselves functions of the regionally varying pro- 
cess coefficients, implying regional specificity is pres- 
ent at a more fundamental level and cannot be 
removed. However, in considering the somewhat 
counter-intuitive features of the optimal weights, as 
described above, it is useful to understand a context 
in which the weights make intuitive sense. 

The specification of the subject-region model can 
be achieved by implementing a stepwise procedure. 
As an initial step, differentials are computed as the 
differences between the subject-region coefficients 
and their respective complement-region weighted- 
average estimates. Each differential is then divided 
by its respective standard errors (see the Appendix 
for a derivation of the covariance matrix for differen- 
tials) to obtain a differential ^-statistic, the distribu- 
tion of which is known to be standard normal in 
large samples under the assumptions that each 
region's residuals are independent and identically 
distributed and that the true differential is zero. If 
the p-value of the differential having the smallest 
absolute value ^-statistic is less than a predetermined 
significance threshold (say, 0.05), then all differen- 
tials are presumed to be statistically significant and 
the stepwise procedure is terminated with the result 
that none of the potential constraints associated with 
the region are imposed in the hybrid model. Con- 
versely, if the minimum absolute value ^-statistic is 
not statistically significant, the corresponding process 
coefficient is constrained to equal its weighted- 
average estimate and the subject-region model is 
re-estimated. Note that only the subject-region model 
is re-estimated, the complement-region fixed-effects 
estimates and implied weighted averages are not 
altered. It is not necessary to revise the complement- 
region fixed-effects estimates because these estimates 
are statistically consistent, meaning that they 
approach their true values as sample size becomes 
large, regardless of the cross-region constraints pres- 
ent in the regional fixed-effects model. Revised differ- 
entials are computed by differencing the re-estimated 
subject-region coefficients from their respective com- 
plement-region weighted-average estimates, the same 
weighted averages that were used in the previous 
step. Standard errors of the revised differentials are 
derived using the methods developed by Murphy and 
Topel (1985) (see the Appendix and Supporting Infor- 
mation for a derivation of their method as it pertains 
to the present application) to account for uncertainty 
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in the value of the constrained coefficient, the coeffi- 
cient that is set equal to the complement-region 
weighted-average estimate. The revised differentials 
are normalized by their respective standard errors 
and the resulting ^-statistics of the unconstrained 
coefficients are ranked to determine the least statisti- 
cally significant differential. If the least statistically 
significant differential has a p-value below the signifi- 
cance threshold, then all remaining differentials are 
statistically significant and the stepwise procedure is 
terminated without imposing any additional con- 
straints; otherwise, the subject-region coefficient cor- 
responding to the least statistically significant 
differential is set equal to the respective complement- 
region weighted-average estimate and the stepwise 
procedure is recycled through another iteration. 

The only exception to the stepwise procedure 
described above concerns cases in which the comple- 
ment-region weighted-average estimate is inconsistent 
with physical constraints imposed on the SPARROW 
model (all source and aquatic removal coefficients are 
constrained to be nonnegative). Invalid weighted-aver- 
age estimates are possible because nonzero weights 
are applied to all complement-region coefficients, not 
just to like-process coefficients. In such cases, the sub- 
ject-region coefficient is not constrained to equal the 
corresponding (invalid) weighted-average estimate, 
even if the associated differential is not statistically 
significant. 

At the conclusion of the stepwise procedure, the 
water manager has determined a re-specified model 
for the subject region, a model that incorporates com- 
plement-region information into the analysis by con- 
straining specific subject-region coefficients to equal 
complement-region weighted-average estimates. Our 
interest with this model, at least with regard to the 
development of a hybrid model, is not the particular 
subject-region estimates; rather, our interest is in the 
particular constraints that are identified. Each con- 
straint represents an empirically valid linear restric- 
tion on the coefficients of the entire regional fixed- 
effects model, the constants of these constraints being 
the weights associated with the weighted-average 
estimates. 

In form, the constraints derived from the "water 
manager's perspective" are more complicated than 
the relatively simple pair-equivalence constraints 
described above for the most general constraint sys- 
tem. Collectively, however, the two systems of con- 
straints are equivalent in terms of potential 
restrictiveness. As shown in the Appendix, if all coef- 
ficients, across all regions, are re-estimated to con- 
form to the constraints implied by all the weights 
across all regions, the fact that the weights sum to 
one is sufficient to cause all like-process fixed-effects 
coefficients to be equal, making the model equivalent 



to the national model. Thus, the set of all possible 
constraints implied by the weighted averages exactly 
duplicates the most restrictive collection of pair- 
equivalence constraints previously described. 

Duplication of the stepwise procedure over all R 
regions provides a complete set of constraints. How- 
ever, as shown in the Appendix, K of the differen- 
tials, across all regions, are linearly dependent, 
implying K of the hypothesis tests undertaken in the 
collective stepwise procedures are not independent. 
Moreover, a logical inconsistency may exist in the 
results of the stepwise procedure applied to all 
regions. To see this, suppose there is only one process 
coefficient in a model consisting of three regions, 
implying there are three fixed-effects coefficients to 
be estimated and three potential constraints to evalu- 
ate using the stepwise procedure, one constraint for 
each region. If the stepwise procedure identifies the 
constraint is binding in the first region and, in the 
second region, the constraint is not binding, then an 
inconsistency arises if the constraint in the third 
region is determined to be binding. Due to the depen- 
dence across constraints, if two constraints are 
binding then the third must also be binding. Unfortu- 
nately, there is no guarantee that application of the 
stepwise procedure to all regions identifies a consis- 
tent set of constraints. 

In the case of inconsistency across constraints, it is 
arbitrary whether to accept all identified constraints. 
Moreover, with multiple process coefficients it is diffi- 
cult to know whether inconsistency is present. One 
way to guarantee there are no inconsistencies in the 
accepted constraints, which also implies that all 
imposed constraints are independent, is to exclude 
one of the regions from the stepwise procedure. The 
exclusion of a region possibly leads to too few con- 
straints; however, too few constraints may be prefera- 
ble to imposing too many, and is consistent with the 
general objective of minimizing the number of con- 
straints applied in the hybrid model (see the above 
discussion pertaining to the use of optimal weights). 

The perspective adopted above, whereby water 
managers specify individual regional models that 
make use of extra-regional information, provides a 
feasible approach to identifying the constraints used 
in the specification of the hybrid model. Given the 
estimates from the unrestricted regional fixed-effects 
model, each regional manager uses a stepwise proce- 
dure to evaluate at most K possible constraints. 
Moreover, the evaluation of each constraint requires 
the estimation of a relatively simple regional model 
consisting of only K coefficients. Duplication of this 
effort over R - 1 regions provides a complete set of 
constraints, implying the total number of evaluated 
constraints is no more than K x (R - 1). If there are 
many regions (R > 3), this number is considerably 
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less onerous than that for the general pair- 
equivalence constraint set noted above. 

As a final step, the empirically validated con- 
straints described above are imposed in the re- 
estimation of the regional fixed-effects model to pro- 
duce the hybrid model. The assumed properties of 
the SPARROW model residuals imply the covari- 
ances of the coefficient estimates are consistently 
estimated. By a well-known theorem (Amemiya, 
1985; theorem 3.2.6), this implies the constants of 
the linear constraints are also consistently estimated 
and, with large samples, there is no statistical 
advantage to re-estimating them in light of the con- 
straints identified by the stepwise process. The con- 
straints imply specific coefficients in the hybrid 
model are restricted to equal a weighted average of 
hybrid model coefficient estimates for the comple- 
ment regions. Although the weights are identical to 
those produced from the original regional fixed- 
effects model, the weighted averages are not, the 
coefficients to which the weights are applied having 
taken on the new values obtained through the simul- 
taneous estimation of all coefficients in the hybrid 
model. Note, also, that due to the nature of the 
weights (see above and the Appendix), any con- 
straint includes coefficients from all regions of the 
model. This means information from every region is 
used in the estimation of any individual coefficient. 
Like the coefficients of the regional fixed-effects 
model, hybrid model coefficients are consistent. The 
imposition of the constraints limits the regional vari- 
ation in hybrid model coefficients but also reduces 
the standard error of their estimates, the estimation 
of which requires a special procedure (Gallant, 1987; 
see the derivation of the standard errors of 
the hybrid model in the Appendix and Supporting 
Information) to account for the constraints. 

In summary (refer to Figure 2), the hybrid model 
is a linearly constrained version of the regional fixed- 
effects model. The coefficient covariances from the 
regional fixed-effects model are used to compute a set 
of weights, these weights representing the constants 
of potential linear constraints to be applied in the 
analysis. To determine which of the potential con- 
straints are actually imposed, the weights are com- 
bined with the coefficient estimates from the regional 
fixed-effects model to obtain a complement-region, 
weighted-average estimate for each fixed-effects coef- 
ficient in the model. A stepwise procedure is used to 
identify which fixed-effects coefficients are signifi- 
cantly different from their associated complement- 
region weighted average. The specific constraints 
applied in the estimation of the hybrid model 
are those associated with weighted averages that 
are not significantly different from a fixed-effects 
coefficient. 



EVALUATION OF RESULTS OF THE HYBRID 
MODELS FOR TOTAL NITROGEN AND TOTAL 
PHOSPHORUS 



As explained in the Introduction, the precision of 
the estimated coefficients of the regional fixed-effects 
model can be improved, with a corresponding loss in 
model fit, by imposing cross-region constraints on the 
regional fixed-effects coefficients, effectively expand- 
ing the scale of the analysis beyond a strictly regional 
scale. To demonstrate the tradeoff between coefficient 
precision and model precision, the methodology of the 
hybrid model described in the preceding section is 
applied to the regional fixed-effects model estimated 
for the three regions shown in Figure 1. 

The stepwise procedure was used to identify which 
of the fixed-effects coefficients to constrain to equal a 
linear function of the complement-region fixed-effect 
coefficients, the parameters of the linear function 
being the optimal weights computed from the regio- 
nal fixed-effects model. As explained in the previous 
section (and shown in the Appendix), the set of all 
possible constraints is not mutually independent, and 
an independent and consistent set of constraints can 
be obtained by excluding one of the regions from the 
stepwise procedure. Even though the selection of 
which region to exclude is arbitrary, the results of 
the analysis are not invariant to which region is cho- 
sen. The East region has the most observations (see 
Figure 1), making this region the likely source of 
regionally specific coefficients and, consequently, 
most likely the source of the fewest constraints. 
Therefore, the East region was excluded from the 
stepwise procedure. Note that the exclusion of the 
East region from the stepwise procedure does not 
mean that the fixed-effect coefficients of this region 
are excluded from the overall set of constraints 
applied in the hybrid model. As remarked above, each 
applied constraint includes every fixed-effect coeffi- 
cient in the model, including the region excluded 
from the stepwise procedure. The inclusion of the 
East region coefficients in each imposed constraint 
implies that the hybrid model estimates of these coef- 
ficients will differ from their values obtained with the 
regional fixed-effects model. 

To evaluate the validity of the £-tests used in the 
stepwise procedure, the residuals from the regional 
fixed-effects model for the Northwest and Southwest 
regions, the only regions subjected to the stepwise pro- 
cedure, were assessed for independence and homosce- 
dasticity (the statistical property that residuals of the 
model have a common variance). For both regions, for 
both TN and TP, the Moran's / was found to be statisti- 
cally insignificant, the smallest p-value among the 
four statistics being 0.38, a value that is not even 



Journal of the American Water Resources Association 



1161 



JAWRA 



Schwarz, Alexander, Smith, and Preston 



marginally significant. The assumption of homoscedas- 
ticity of the residuals was evaluated using a variant of 
a test for heteroscedasticity developed by White (1980), 
whereby the squared residuals (for TN, the squared 
residuals are multiplied by the observational weights) 
of the regional fixed-effects model are regressed on the 
derivatives of the model predictions at each monitored 
reach with respect to each model coefficient (i.e., the 
gradients of the nonlinear SPARROW model - the 
analog of the explanatory variables in a linear model). 
In all cases, the F-statistics for these regressions 
were not significant at the 0.05 level, indicating that 
heteroscedasticity was not a concern. 

The stepwise procedure was used to identify 22 
constraints, each constraint being applied to the 54 
fixed-effects coefficients (41% of the coefficients were 
constrained) included in the TN model. The remain- 
ing 32 potential constraints were not applied either 
because they pertain to potentially nonindependent 
constraints associated with the East region (18 coeffi- 
cients), were not constrained due to the correspond- 
ing weighted-average estimate being inconsistent 
with physical bounds that all source and aquatic 
removal rate coefficients be nonnegative (2 coeffi- 
cients in the Southwest region), or had statistically 
significant regional differentials (12 coefficients, 7 in 
the Northwest region and 5 in the Southwest region). 
Coefficients that were unconstrained, aside from all 
coefficients in the East region, were: wheat and pas- 
ture/range nitrogen application, barren land, artifi- 
cial drains, and stream removal rate (both the 
Northwest and Southwest regions), corn/soybean 
nitrogen application and temperature (Northwest 
region), and urban population and soil permeability 
(Southwest region). For the TP model, the stepwise 
procedure identified 28 constraints to be applied to 
the 45 fixed-effects coefficients (62%), consisting of all 
coefficients for both the Northwest and Southwest 
regions, except the Southwest region values of the 
specific catchment area and stream removal rate 
coefficients, both of which had statistically significant 
differentials. 

The TN and TP regional fixed-effects models were 
re-estimated with nonnegative physical bounds placed 
on all source and aquatic removal rate coefficients, 
and with the 22 and 28 cross-region linear constraints 
identified from the stepwise process. Statistical sum- 
maries of the resulting estimations of the hybrid mod- 
els, and their comparison to statistics for the national 
and regional fixed-effects model without cross-region 
constraints, are given in Table 2. In terms of percent 
of coefficients that are statistically significant at the 
0.05 level, the hybrid models perform much better 
than the regional fixed-effects models, for both TN 
and TP - a consequence of the cross-region con- 
straints, which greatly increase the error degrees of 



freedom (number of observations, plus number of con- 
straints, minus number of estimated coefficients) of 
the hybrid model as compared to the regional model. 
The TN hybrid model average coefficient of variation 
(the ratio of each coefficient's standard error to its 
estimated absolute value, averaged over all coeffi- 
cients in the model), a normalized measure of the 
imprecision of all the coefficients in the model, is 
about 38% [100 (0.891 - 0.549)/0.891] lower than that 
for the regional fixed-effects model. The hybrid model 
improvement in coefficient precision is even greater 
for TP, with an average coefficient of variation that is 
84% [100 (2.416 - 0.389)/2.416] lower than that for 
the regional fixed-effects model. 

The reduced model fit caused by the imposition of 
constraints to obtain the hybrid model can be mea- 
sured by a number of statistics, including the RMSE, 
R 2 , and yield R 2 . As the model is estimated in loga- 
rithmic space, 100 times the RMSE can be inter- 
preted as the approximate standard error of a 
predicted load for any arbitrary reach in the model, 
expressed as a percent of the predicted load (this 
approximation is derived from a first-order Taylor 
expansion of the error in load, expressed in real 
space; see Schwarz etal., 2006, for a related discus- 
sion). The R 2 and yield R 2 are computed from load 
and predicted load having a logarithm transforma- 
tion. The yield R 2 - the R 2 computed using predicted 
and observed values of the logarithm of load divided 
by upstream basin area - is a measure of fit that 
effectively removes from R 2 the inflating effects from 
estimating a model that encompasses a wide range of 
basin scales (Schwarz et al., 2006). 

In terms of RMSE, the TN hybrid model is only 
slightly less precise than the regional fixed-effects 
model, the difference in RMSE between the hybrid 
and regional fixed-effects models being only 0.007 
(0.505 - 0.498). Thus, using the approximation des- 
cribed above relating differences in RMSE to differ- 
ences in the standard error of load, expressed as a 
percent of load, the predicted TN load for the hybrid 
model has a standard error that is only 0.7% greater 
than that of the regional fixed-effects model and 
nearly 4.8% [100 (0.505 - 0.553)] lower than the 
national model. Similarly, the R 2 and yield R are 
nearly equal between the regional fixed-effects and 
hybrid TN models, the difference being only 0.008 for 
R 2 and 0.015 for yield R 2 . 

The TP hybrid model has relatively more con- 
straints applied than the TN model, implying the sac- 
rifice in model fit as compared to the regional fixed- 
effects model is greater. The RMSE of the hybrid 
model is approximately halfway between the national 
model and the regional fixed-effects model, with the 
standard error of predicted load, expressed as a 
percent of load, being approximately 3.7% [100 
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(0.726 - 0.689)] greater than that for the regional 
fixed-effects model and 3.6% [100 (0.726 - 0.762)] less 
than that for the national model. The R 2 and yield R 2 
of the TP hybrid model are only 0.019 and 0.046 less 
than those for the regional fixed-effects model. 

The relative model fits of the hybrid and regional 
fixed-effects models can also be evaluated in terms of 
an F-test applied to the sum of squared errors and 
error degrees of freedom statistics reported in 
Table 2. For TN, the p-value of the F-statistic evalu- 
ating the statistical significance of the difference in 
sum of squared errors between the two models is 
0.061, indicating that the model fit of the regional 
fixed-effects model is not significantly better than the 
hybrid model. The F-test for the TP model gives a 
p-value <0.001, implying the regional fixed-effects 
model has a significantly better fit than the hybrid 
model. 

Perhaps a more useful measure of model perfor- 
mance is prediction accuracy. Such an assessment 
includes not just model error (as expressed by the 
MSE), but also the sample error associated with 
uncertainty in the coefficient estimates. To assess 
model accuracy we compiled concentration and flow 
data (see Saad et al., this issue) to compute mean 
loads, detrended to a base year of 1992 (see the load 
estimation method described in Alexander et al. , 
2008), for an additional 678 (TN) and 865 (TP) sta- 
tions linked to the E2RF1 reach network but not used 
in the estimation of any of the models described in 
this article. Predictions from the national, regional 
fixed-effects, and hybrid models were used to compute 
error as a percent of the monitored mean load values 
at the reach locations of the additional stations (see 
Schwarz et al., 2006, for a description of the predic- 
tion methodology). Table 2 reports the median abso- 
lute value of the percent prediction error for TN and 
TP across all stations for each model. The results are 
very similar to those obtained for RMSE. For TN, the 
three models have nearly the same accuracy, with 
the national model being least accurate (median abso- 
lute prediction error of 42.3%) and the hybrid model 
having slightly better accuracy (38.4% absolute error) 
than the regional fixed-effects model (40.2% absolute 
error). Predictions from the TP models are less accu- 
rate than the TN models, with the national TP model 
having the least accuracy (67.1% absolute error) and 
the regional fixed-effects model the greatest accuracy 
(52.6% absolute error). As with RMSE, the TP hybrid 
model has a median absolute error (59.1%) that is 
approximately halfway between the national and 
regional fixed-effects models. 

The Moran's / statistic, derived from the residuals 
of the hybrid model and reported in Table 2, is statis- 
tically significant for both TN and TP, indicating the 
presence of spatial correlation in these models. Given 



the detection of spatial correlation in the regional 
fixed-effects model, this result is not unexpected for 
TP. The presence of spatial correlation for the TN 
hybrid model, however, suggests that the constraints 
applied to obtain the TN hybrid model have induced 
a regional pattern in the residuals. 

The coefficient estimates, standard errors and 
^-statistics for the TN and TP hybrid models, along 
with detailed results for the national and regional 
fixed-effects models are reported in the Supporting 
Information. Figures 3 and 4 summarize these results 
in a way that allows for a visual comparison of the 
estimates. The figures display the coefficient esti- 
mates for the hybrid and regional fixed-effects models 
expressed as a percent difference from the national 
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FIGURE 3. A Comparison of Regional Fixed-Effects and Hybrid 
Model Estimates for Total Nitrogen (TN). Models are based on data 
from Alexander et al. (2008). Plotted values represent the coeffi- 
cient estimate expressed as a percent difference from the national 
model estimate (see Alexander et al., 2008, for the national model 
coefficient estimates). Circled values pertain to regional fixed-effect 
and hybrid model coefficients that are constrained in the stepwise 
procedure and consequently generate a constraint in the estimation 
of the hybrid model. Note that the horizontal scale is compressed 
in both the right and left margins, so much so that values on the 
left and rightmost reference lines are unresolved. 
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FIGURE 4. A Comparison of Regional Fixed-Effects and Hybrid 
Model Estimates for Total Phosphorus (TP). Models are based on 
data from Alexander et al. (2008). Plotted values represent the coef- 
ficient estimate expressed as a percent difference from the national 
model estimate (see Alexander et al. , 2008, for the national model 
coefficient estimates). Circled values pertain to regional fixed-effect 
and hybrid model coefficients that are constrained in the stepwise 
procedure and consequently generate a constraint in the estimation 
of the hybrid model. Note that the horizontal scale is compressed 
in both the right and left margins, so much so that values on the 
left and rightmost reference lines are unresolved. 

model estimates. Thus, a value of zero on the graph 
implies the coefficient equals the national model esti- 
mate, and hybrid and regional model estimates that 
have a binding physical constraint (zero values for 
source and aquatic removal rate coefficients con- 
strained to be nonnegative) are plotted on the -100% 
difference reference line (physical constraints imply 
no hybrid or regional model difference estimates are 
less than -100% for any source or aquatic removal 
process coefficient). Although not readily apparent in 
all cases, each process coefficient includes three plot- 
ted values, one value for each region, for each of the 
hybrid and regional fixed-effects models. Circled val- 
ues pertain to hybrid and regional model coefficients 
that are identified in the stepwise procedure to be 
constrained to equal a weighted average of the com- 
plement-region coefficient estimates; thus, hybrid 



model values that are not circled are unconstrained. 
In order to display all the coefficient estimates in a 
single graph it is necessary to use a nonlinear scale, 
one that greatly compresses the horizontal scale in 
the positive and negative tails of the figure, so much 
so that values less than -800% or >800%, plotted 
along the left and rightmost reference lines, are effec- 
tively unresolved. For TN, values along the rightmost 
reference line range from 1,499% (hybrid model artifi- 
cial drains estimate in the Southwest region) to 
4,760% (hybrid model barren estimate in the North- 
west region); for TP, values near or along the right- 
most reference line range from 773% (regional fixed- 
effects model barren estimate in the Southwest 
region) to 2,402% (regional fixed-effects model shrub 
estimate in the Northwest region). 

A comparison of the various model coefficient esti- 
mates for TN is shown in Figure 3. Generally, the 
range of regional variability differs considerably across 
the process coefficients, with the national model pro- 
viding a representative estimate for each coefficient. 
As expected, the greatest range of variation, for most 
process coefficients, is exhibited by the regional fixed- 
effects model estimates. Large (absolute value >200%) 
differences from the national estimates are observed 
for the regional fixed-effects estimates of a number of 
process coefficients, including: barren and shrub lands, 
artificial drains, soil permeability, temperature, spe- 
cific catchment area, drainage density, and reservoir 
removal. Variation in the hybrid model estimates 
approximates the variation in the regional fixed-effects 
estimates for a few process coefficients (wheat, pas- 
ture/range, barren and shrub lands, artificial drains, 
soil permeability, temperature, and stream removal); 
generally, however, the hybrid model estimates are 
closer to the national model estimates than are the 
regional fixed-effects estimates. Two hybrid model 
coefficients that represent the largest differences from 
the national model estimate, shrub land (East region) 
and barren land (Northwest region), correspond to the 
smallest sources in their respective regions (Alexander 
et al., 2008), so the effect of the large coefficients on 
predicted load is relatively minor. Using ±100% differ- 
ence from the national model estimate as a benchmark 
for a substantial regional variation, there are 10 such 
coefficients (19% of the 54 coefficients) exceeding this 
threshold for the TN hybrid model. Hybrid model esti- 
mates subject to a cross-region constraint (i.e., hybrid 
estimates that are constrained to equal a weighted 
average of the complement regions' hybrid model esti- 
mates; values that are circled in the figure) are gener- 
ally within 100% of the national model estimates, 
except for the drainage density and reservoir removal 
estimates. 

A comparison of coefficient estimates between the 
national, regional fixed-effects, and hybrid TP models 
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is shown in Figure 4. As with TN, the national model 
provides a representative estimate for each coefficient. 
Unlike TN, the TP regional fixed-effects model esti- 
mates show greater variation than the hybrid model 
estimates for all process coefficients. The variation in 
the hybrid model estimates approximates that of the 
regional fixed-effects estimates only for shrub land and 
stream removal, the latter case being an example 
where both the regional fixed-effects and hybrid model 
estimates face a binding physical constraint (both esti- 
mates are at plotting position -100%, implying they 
are evaluated at zero - the physical bound placed on 
aquatic removal rate coefficients). Many of the most 
extreme regional fixed-effects estimates in Figure 4 
appear spurious, their values being constrained in the 
hybrid model estimation (circled values). Only the coef- 
ficient for shrub land has a large deviation from the 
national model and is not subsequently constrained in 
the hybrid model estimation, although this value is 
unconstrained because it pertains to the East region, 
the region excluded from the stepwise process. As 
remarked above, shrub land pertains to an insignifi- 
cant source in the East, so the effect of the coefficient 
on predicted load is negligible. All hybrid model esti- 
mates subject to a cross-region constraint (the circled 
values) lie within ±100% of the national model esti- 
mate, which implies that many of the extreme regional 
fixed-effects estimates become "controlled" by the 
hybrid model methodology. If the benchmark for a sig- 
nificant regional variation is an absolute difference 
greater than 100% from the national model estimate, 
the TP hybrid model has only one such coefficient (the 
shrub land coefficient for the East region; representing 
only 2% of the 45 total coefficients). 

The results of the analysis show the hybrid 
approach achieves greater precision in the coefficient 
estimates as compared to the regional fixed-effects 
model, with only a modest sacrifice in prediction 
accuracy. The greater precision in the hybrid model 
coefficients for TN comes at virtually no cost in terms 
of model fit, whereas the much larger improvement 
in the precision of the TP hybrid model coefficients 
requires a sacrifice of approximately half of the model 
fit advantage of the regional fixed-effects model over 
the national model. Very few hybrid model coeffi- 
cients show evidence of large regional variation, the 
most notable cases involving barren and shrub land 
in both the TN and TP models. It is possible that the 
data for these land classes contain classification 
errors, or the contaminant sources they represent are 
not sufficiently differentiated, and an improved model 
specification would be possible by combining the two 
classes into a single variable. Both the TN and TP 
hybrid models exhibit spatial correlation in the resid- 
uals, a likely byproduct of the coarse regionalization 
adopted for this study. 



DISCUSSION 



An issue frequently raised by the regionalization 
approach concerns the source of regional variation in 
the coefficients. From a modeling perspective, it is 
reasonable to suppose that the causes of regional var- 
iation are unrelated to variables that could be 
employed in the statistical analysis, for if such vari- 
ables were the cause their use in the model would 
alleviate the need for regionalization. Consequently, 
regionalization likely reflects the absence of relevant 
data in the analysis. For example, the natural phos- 
phate content of basin minerals, a prominent variable 
affecting TP (Garcia et al., this issue), is absent from 
the models described here due to a lack of consistent 
data at the national scale. 

Even if pertinent variables are available, regional 
patterns in the residuals may be evident if the vari- 
ables are measured with error, or if the functional 
relation between water quality and the causative 
variables is too complex to be adequately specified in 
a large-scale model. An example of a variable mea- 
sured with error from the TN model is atmospheric 
deposition, which pertains only to wet deposition and 
excludes dry deposition. In addition, the surrogate 
land use variables employed in both the TN and TP 
models are incomplete measures of nutrient loading 
for specific land types, and are likely the source of 
regional coefficient variation. Given the scale of the 
analysis, therefore, regional variability in the coeffi- 
cients is expected to be pervasive, making it difficult 
to attribute a specific cause. 

The fixed-effects approach is just one of many that 
have been used to introduce regional variation into 
statistical models. An alternative is the "random 
effects" approach (Judge et al., 1985), also called the 
hierarchical or multilevel method (Reckhow et al. , 
2009; Kashuba et al, 2010; Qian et al., 2010), which 
specifies model coefficients to be random variables 
having an assumed multivariate probability distribu- 
tion with unknown parameters. The advantage of 
this approach is that only the parameters of the mul- 
tivariate distribution for the coefficients need to be 
estimated, not the individual effects for each region, 
the estimation of which could be difficult if there is 
insufficient within-region variation in basin attri- 
butes. The multilevel method also has the potential 
to provide a richer description of the origin of error 
in the model, which is useful in the evaluation of pre- 
diction uncertainty. However, the method places 
rather strict assumptions on the statistical distribu- 
tion of the regional coefficients (Mundlak and Yahav, 
1981); typically, the regional coefficients are assumed 
to be normally or log-normally distributed and inde- 
pendent across regions. Moreover, the method may be 
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difficult to implement computationally in the context 
of a nonlinear model, such as SPARROW, particu- 
larly if the multivariate distribution of the coeffi- 
cients has a complicated covariance structure. 

Another popular approach is the region of influ- 
ence approach (Acreman and Wiltshire, 1987; Burn, 
1990a,b), which posits a different set of coefficients 
for each observation in the analysis. The procedure is 
to develop a unique estimation sample for each obser- 
vation, one that has explanatory variable characteris- 
tics that come closest to the given observation. 
Similar methods are the regression tree analysis and 
the cluster analysis (Burn and Boorman, 1993; 
Robertson et al., 2006; Soranno et al., 2010). A com- 
mon advantage of these "data-driven" methods is the 
flexible concept used to define a region: the delinea- 
tion of a region is based on all model variables, not 
solely on the spatial dimension. However, this flexi- 
bility also makes it difficult to interpret the estimated 
coefficients and compare them to values obtained in 
other studies. For this reason, these methods may be 
better suited for prediction of water-quality charac- 
teristics at unmonitored sites; they may be less 
appropriate for applications requiring precise esti- 
mates of individual model coefficients, as is necessary 
for the evaluation of water-quality response to alter- 
native management or forecasting scenarios. 

A common purpose for employing regionalization 
methods is to account for regional variation in the 
data even where the causes of the variation are diffi- 
cult to identify explicitly. The fixed-effects method, 
like the other methods described in this section, 
assigns regional variation to specific process coeffi- 
cients of the model, an assignment that may provide 
some insight as to the underlying latent biogeochemi- 
cal or hydrological processes, human activities, or 
data imperfections driving regional variation (see 
Hoos and McMahon, 2009, for a discussion of how a 
causal relationship could be identified in the context 
of a fixed-effects approach to regionalization). How- 
ever, the lack of a complete causal understanding for 
regional variation with these approaches may limit 
the utility of regionalized models in the generation of 
model predictions under alternative management sce- 
narios. For example, if regional variation reflects dif- 
ferences in unspecified climate conditions, the 
regional model may be less useful for simulating the 
effects on water quality resulting from climate 
change - without knowledge that the regional effect 
is due to climate it is not possible to modify this effect 
for different climate conditions. This issue also has 
implications for quantifying the uncertainty of simu- 
lations based on the model. The inclusion of regional 
effects may improve model fit, but given the lack of 
understanding regarding the nature of the regional 
variation, it remains an open question as to whether 



this reduced error should translate into lower error 
in the simulation results. 

The fixed-effects approach leaves unanswered a 
clear and definitive methodology for defining the geo- 
graphic limits of a region. The methodology imposes 
a generally arbitrary regional structure, usually 
defined according to hydrological properties (e.g., 
watersheds), and then uses statistical evidence to 
determine the relevance of this structure. Other 
methods, such as cluster and regression tree analyses 
and the method of regional regression, offer a more 
flexible approach that uses characteristics of the data 
to define regions. However, if the underlying factors 
causing regional variation are unknown, uncertain- 
ties will remain with these approaches as to how 
reliably to use regional models in simulation applica- 
tions. Another approach may be to delineate regions 
based on a preliminary assessment of residuals using 
a nonregionalized model. As is shown in Figure 1, 
residual patterns are evident for both TN and TP 
along the coastal areas of the eastern U.S. and Gulf 
of Mexico, patterns that do not fit perfectly with the 
regional basins used for this study. However, as 
remarked above, each coefficient in a SPARROW 
model must interact with a spatially varying land or 
aquatic causal variable, so that regional patterns in 
the model coefficients need not translate into discern- 
able regional patterns in the residuals, drawing into 
question the utility of this approach. 

The standard approach to regionalization is one in 
which each regional model is estimated indepen- 
dently (as with the NAWQA MRB regional models 
included in this Featured Collection). As discussed in 
the Introduction, the principal disadvantage of this 
approach is insufficient within-region spatial varia- 
tion of the explanatory variables, leading to imprecise 
estimates of the model coefficients. There are, how- 
ever, several mitigating circumstances acting in favor 
of such an approach. The independent regional analy- 
sis may allow for scaling the level of effort required 
to develop certain geospatial datasets. For example, 
the resources required to produce a good point source 
dataset should scale, at least approximately, with the 
basin area of the analysis. It may also be possible in 
a regional analysis to focus resources on the develop- 
ment of certain geospatial datasets that are impor- 
tant for the region of interest but of lesser 
importance elsewhere. A map showing the location of 
millponds might be of particular interest for a model 
pertaining to the Northeast U.S., but of lesser impor- 
tance for a model focused on the Midwest. The 
independent regional analysis also allows greater 
flexibility in the specification of the model, such as 
the capability to use geospatial data that are avail- 
able only for the region of interest. This latter point 
brings to attention one of the weaknesses of the 
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fixed-effects approach, that it does not address differ- 
ent specifications of the model, either in terms of 
functional forms assumed for the transport processes 
or the particular contaminant sources and other vari- 
ables used in the analysis. Different specifications of 
the model across different regions may be required to 
address cases in which a source is poorly defined in a 
specific region leading to a potentially spurious esti- 
mate of a model coefficient, such as the cases noted 
above in the results of the hybrid models for the bar- 
ren land TN model coefficient in the Northwest 
region and the TN and TP model shrub land coeffi- 
cients in the East region. The independent NAWQA 
MRB models of this Featured Collection demonstrate 
the utility of incorporating into their analyses factors 
other than those used in a national analysis. Indeed, 
differences in these factors, as well as a different 
time period (the independent regional models use a 
base year of 2002, compared to the year 1992 used 
for this study) and a large difference in the number 
of stations used in model estimation may help explain 
any discrepancies in model results among this study 
and the independent MRB models. 

The results of this study show that there are likely 
to be gains in model inference by combining the results 
of independent regional models with those of national- 
scale models, a direct consequence of more information 
being better than less. Perhaps a better incarnation of 
this maxim, however, is achieved through Bayesian 
analysis. Under a Bayesian approach, the national- 
scale model could serve as prior information for the 
estimation of a regional scale model (Qian et al, 2005; 
Qian and Reckhow, 2007). Rather than deciding 
between coefficients based on information from either 
the national or regional scale, as is done under the 
fixed-effects approach described in this article, the 
Bayesian posterior distribution for the regional coeffi- 
cients would reflect the combined information from 
both scales. 

The method of regionalization employed in this 
article, whereby the variation of coefficients across 
space is constrained to improve statistical inference, 
can also be applied to the analysis of the variation of 
coefficients across time. In a temporal context, one 
can imagine multiple SPARROW models, each cover- 
ing the same study area, estimated independently 
according to water-quality conditions at different 
points in time. The pertinent question is whether it is 
necessary to specify different values for the coeffi- 
cients for each modeled time period, or whether a 
given coefficient could be constrained to equal a 
weighted average of the coefficients for the other peri- 
ods. Conceivably, the method described here could be 
used to form cross-period constraints on the coeffi- 
cient estimates. The method may need to be modified 
if the assumption that the errors in the coefficient 



estimates are independent across time is untenable, 
as would be the case if the residuals of the 
SPARROW model were serially correlated. 



SUMMARY AND CONCLUSIONS 



The findings of this study demonstrate that region- 
alization of national SPARROW models for TN and 
TP using the fixed-effects approach has the potential 
to provide modest, albeit highly statistically signifi- 
cant improvements of 6 and 7%, respectively, in the 
accuracy of model predictions. However, this 
enhanced accuracy comes at a cost of reduced preci- 
sion of the regional fixed-effects coefficient estimates 
relative to the national model estimates. A method is 
proposed to combine national and regional informa- 
tion to expand the variation in the underlying process 
variables, thereby improving the precision of the 
regional fixed-effect coefficient estimates. The method 
is applied to three regions delineating major water 
basins of the conterminous U.S., leading to the speci- 
fication of hybrid models for TN and TP. The results 
of these analyses show that the method is successful 
in improving the precision of the process coefficient 
estimates as compared to a standard regionalized 
fixed-effect approach, while providing modest reduc- 
tions in the standard error of TN and TP model load 
predictions, as compared to the national model with- 
out regionalization, of 5 and 4%. The hybrid model 
estimates show little evidence of strong regional vari- 
ation in the coefficient estimates. Only 19% of the TN 
hybrid model coefficients and just 2% of the TP 
hybrid model coefficients have more than ±100% devi- 
ation from the national model estimate. This lack of 
strong regional variation is possibly due to the coarse 
partitioning of the U.S. into just three regions. 

From a theoretical perspective, the regional bound- 
aries used in a regionalization analysis should delin- 
eate natural areas having approximately homogeneous 
values of the underlying model coefficients. As the 
hybrid model developed for this study constrains the 
regional variation of individual coefficients, all coeffi- 
cients need not adhere to identical spatial patterns of 
homogeneity. This also means the method can be use- 
ful in limiting the variation of coefficients across 
regions that do not strictly conform to the natural 
areas described above, such as the case where the 
regional analysis requires the delineation of areas 
along arbitrary political boundaries to accommodate 
specific management considerations. 

The national model residuals for TN and TP display 
distinctive regional patterns; however, these residual 
patterns may not be fully indicative of spatial patterns 
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in the underlying values of the model coefficients. 
Additional work is required to better define regions of 
the nation having the characteristics of natural areas 
reflecting homogeneous values for model coefficients. 
Such an effort will be greatly aided by including the 
large set of monitoring data compiled as part of the 
NAWQA MRB studies appearing in this Feature Col- 
lection. Those data should allow for the specification of 
many more regions, resulting in greater improvements 
in model fit than obtained here. The added observa- 
tions would also reduce the standard errors of the 
regional fixed-effects coefficients, or at least moderate 
the increase in standard errors caused by the larger 
number of fixed effects associated with a finer regional 
delineation of the nation. Subsequent analysis should 
also investigate alternative formulations of the model 
beyond that of the previously published national mod- 
els to determine if some of the regional patterns in the 
hybrid model coefficient estimates documented above 
are an artifact of model misspecification. Regardless of 
the specific modifications to the analysis enabled by a 
larger dataset, the pooling of information across 
regions, using methods akin to those described in this 
article, should prove helpful in achieving a favorable 
tradeoff between precision of the coefficient estimates 
and precision of the model predictions. 



APPENDIX 



The following presents a mathematical formulation 
of the hybrid model, beginning with a derivation of 
the complement-region weighted-average estimates 
and region-specific coefficient differentials. Let there 
be R regions: the subject region, s, and the R - 1 
complementary regions, indexed by r. Define f5 r to be 
the estimated vector of K process coefficients associ- 
ated with region r. The estimated region s coefficient 
differential, 5 S , is defined as the difference between 
the region s estimated coefficient vector, f) s , and a 
weighted average of the estimated coefficient vectors 
for the complementary regions, 

ds = K-^yv s X (Ai) 

where W s>r is a K x K matrix of weights, with 
J2r^s Ws,r = Ik, Ik being the if-dimensional identity 
matrix. Under the assumption that the error terms of 
the regional models are not spatially correlated 
across regions, the regional coefficient estimates are 
independent, having region-specific covariance matrix 
V r . Thus, the covariance matrix for the region s 
coefficient differential vector is 



V(S S ) =Vs + J2 W ^- V ' W v- ( A2 ) 

Proposition 1: the covariance matrix of the esti- 
mated regional differentials is minimized (i.e., the 
quadratic form x'V^<5 s ^x is minimized for arbitrary 
if x 1 vector x) if the weights are set according to 

Proof: Consider setting the weights equal to the 
proposed value, plus some arbitrary matrix M s>r , so 
that 

W s/ = V/'l V r 1 + M s,r- (A4) 

Because the weights must sum to \ K , it must 
be the case that J2r^s M s ,r = Ok, 0 k being a K x K 
zero matrix. Upon substituting Equation (A4) into 
Equation (A2), the covariance matrix of the regional 
differentials becomes 

+ ( E M '°r) ( E v ' :1 ) + E M-V^M',,, 

= V s + (j2 V ' Tl ) + E Ms -r y r ' M ^- 

(A5) 

As V r is positive definite, M sr V^ 1 M' sr must be a 
nonnegative matrix. Thus, any quadratic form derived 
from Equation (A5) is minimized if and only if 
T, r ^s M s,rVr lM 's,r = Ok- By theorem 12.2.5 of Graybill 
(1969), the sum of nonnegative matrices is zero if and 
only if each of the nonnegative matrices comprising 
the sum are also zero, implying M s r V ; T 1 M' s . r = 0# for 
all r. As V r is positive definite, there exists a nonsingu- 
lar matrix P r such that V; 1 = P r P' r (Graybill, 1969, 
theorem 12.2.2), implying M Sir V^ 1 M' Sir = U s , r U' s , r , 
where U Sir = M sr P r . The ith diagonal element of 
U s ,rU' s .,- is given by J2k=i k> which equals zero for 
each i element if and only if U Str i t k = 0 for all i, k. 
Thus, U s>r = M SjF P r = 0 K , which, given that P r is non- 
singular, implies M s r = 0^. Therefore, Equation (A3) 
is the unique value for the weights that minimizes 
any quadratic form based on the covariance matrix of 
the regional coefficient differentials. 

With weights set according to Equation (A3), the 
covariance matrix for <5 S , is given by 
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(A6) 



An implication of Equation (A6) being a minimum 
covariance matrix is that the diagonal elements have 
the smallest possible value subject to the constraint 
that the weights sum to one. Consequently, the test 
for significance of a given regional differential is 
made as sensitive as possible. Note, also, that the 
assumption of independent errors across regions 
implies that the weights minimize both the variances 
of the differentials and the variances of the comple- 
ment-region weighted-average estimates. 

An examination of some special cases can enhance 
understanding of the properties of the optimal 
weights. The optimal weights described above have a 
straightforward large-sample interpretation if it is 
assumed that the asymptotic covariance matrix of 

\/Nr(J} r — /? r j is identical for each region, where N r is 
the number of monitoring sites in region r. In this case, 
lim VrVj 1 =Ji m {Nj/N r )lK and the asymptotic opti- 



2V 



mal weights are simply the share of total observations 



(5 



in each complement region, W sr 

(N r /E Ms Nj)lK. 

The optimal weights also have a straightforward 
interpretation if the underlying model is linear. Let 
the water-quality model for region r be given by 
y r = X r /} r + e r , where y r is an N r x 1 vector of water- 
quality load measurements, X r is an N r x K matrix of 
explanatory variables, and e r is an N r x 1 vector of 
independent, identically distributed residuals having 
mean zero and variance of. In this case, the differen- 
tial is computed as 

(A7) 

which is the difference between the ordinary least 
squares estimate of the process coefficients for the 
subject region s and the weighted linear least squares 
estimate of the process coefficients for the combined 
R - 1 complementary regions. Other than accounting 
for differences in the mean squared error (MSE) 
across regions, the optimal weights result in a test 
statistic that does not distinguish between individual 
complementary regions and is equivalent to compar- 
ing the ordinary least squares estimate of region r to 
a single like-process coefficient estimated from all 
complement region data pooled together. 

The analysis is more complicated if the underlying 
model is nonlinear in the coefficients, such as in a 
SPARROW model. In that case, the regional covari- 
ance matrices, V r , are functions of the underlying 



regional coefficients, /? r , implying it is not possible to 
simply pool all complement region data. Nevertheless, 
the regional deviation statistic, and its associated 
covariance matrix given in Equation (A6), provides a 
useful ^-statistic for evaluating the significance of 
regional differences in process coefficients. The 
nonlinear estimate of <5 S is asymptotically efficient (in 
the sense of having minimum variance), and the 
^-statistic, computed in conjunction with the associ- 
ated covariance matrix, itself formed from standard 
nonlinear coefficient covariances V, for the individual 
regions (see, for example, Amemiya, 1985), is asymp- 
totically distributed standard normal as long as the 
underlying regional coefficients have the same 
asymptotic properties and sampling of the regions is 
proportionate to total sample size. 

The stepwise procedure applied to subject region s 
requires repeated estimation of the region s model, at 
each step evaluating the statistical significance of the 
regional differentials. The initial evaluation of the 
differentials uses the standard errors implied by 
Equation (A6). Evaluations in subsequent steps per- 
tain to a partially constrained subject-region model in 
which some of the coefficients are set equal to a com- 
plement-region weighted-average estimate. It is not 
necessary to re-estimate the weighted averages at 
each step; the regional fixed-effects model produces 
statistically consistent estimates of the coefficients 
and associated covariance matrices, regardless of any 
underlying cross-region constraints on the coeffi- 
cients. The only complication in the procedure con- 
cerns the standard errors of the unconstrained 
coefficients in the constrained subject-region model, 
which depend on the uncertainty in the values of the 
constrained coefficients (Murphy and Topel, 1985) 
caused by uncertainty in the weighted averages (see 
the inverse term in the right-hand side of Equation 
A6). Using results from Murphy and Topel, it can be 
shown (see the derivation contained in the Support- 
ing Information) that the covariance matrix for the 
vector of regional differentials for subject region s 
corresponding to the unconstrained coefficients of 
that region (denoted 5f) is consistently estimated by 



V 5" =V 



. u 



K °i ) y s V S [U,C\ V 



(b.) 
(b.) 
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(A8) 
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where V s is the weighted nonlinear least squares 
(WNLS) estimate of the covariance matrix for the 
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unrestricted coefficients, the constrained coefficients 
being set equal to the corresponding value of the 
weighted average of the complement-region coeffi- 
cients, denoted here as b s ; of is the MSE of the region-s 
model prior to imposing any constraints; df 1 is the 

MSE of the region-s model with some of the coefficients 

* u 

constrained (the model used to generate V s ); l K u is a 
K.u x K.u identity matrix, being the numbe^ of 
unconstrained coefficients in the region-s model; V s \u,c\ 
is 3lK.u x K% submatrix of the inverse of the covariance 
matrix of the region-s model prior to imposing any 
constraints, the rows of which correspond to the 
unrestricted coefficients and the = K — col- 
umns correspond to the coefficients that have been 
constrained to equal the corresponding value of the 
weighted average of the complement-region coeffi- 



cients; and V 



(6.) 



[AM 



is the submatrix of the estimated 



covariance matrix for the weighted average of the com- 
plement-region coefficients for region s (i.e., the sub- 
matrix of the estimate of (j2j^ s Vjj > see Equation 

(A6) above), the rows and columns of which pertain to 
the group of coefficients identified respectively by A 
and B, where A and B equal either U, for the uncon- 
strained coefficients of the region-s model, or C, repre- 
senting the constrained region-s coefficients. Note that 
the regional weights selected to minimize the covari- 

IX 



ance matrix V(bsJ imply the quadratic form x'V ( <5 
is also minimized. 

The collection of cross-region constraints, identified 
by the application of the stepwise procedure for each of 
R - 1 regions comprising a national model, can be 
applied toward the simultaneous estimation of the 
regional fixed-effects model, producing the hybrid 
model. Mathematically, the collection of constraints 
can be described as follows. Let <5 = {8[, . . . , 8' R } repre- 
sent the vector of stacked regional differential vectors 
for each region, let ji = {f>[, . . . , 0 R } represent the vec- 
tor of stacked regional fixed-effects coefficient values, 
and define G to be a RK x RK matrix that relates the 
stacked regional differentials to the regional fixed- 
effect coefficients, such that 8 = Gfi. From the defini- 
tion of (5, given in Equation (Al), this implies 



Ik -W 1i2 
-W 2 ,i 1 K 



-w ljB 
-w 2jB 

Ik 



(A9) 



where the W S;7 . are defined in Equation (A3). Note 
that the restriction that the weights sum to one implies 
the rightmost column of block matrices in Equation 
(A9) equals minus one times the sum of the R - 1 
leftmost columns, so G does not have full rank and is 



not invertible. A full rank formulation of the first 
R - 1 sets of regional differentials 8 = {8[, . . . , 8' R -i}' 
can be obtained by modifying the /} vector to consist 
of R - 1 sets of coefficient differences from the region 
R coefficient vector, so that 8 = GAjS, where 
A/? = {j8£ -/%,... , p' R _ x - ji' R }' and G consists of the 
(R - 1) x (R - 1) upper-leftmost blocks of G, 



Ia 

-w 2il 



1.1 



-Wl, 2 

Ik 



-Wi^_i 
-W 2ifl _i 



(A10) 



The .Rth set of regional differentials can also be 
expressed as a linear combination of A/i, such that, 
8 R = g R A/J, where g R = [ - W fl ,i . . W flJ? _i ]. As G is 
nonsingular, we have 8 R = g R G 8, and the i?th set of 
regional differentials are shown to be a linear combi- 
nation of the other R - 1 sets of differentials. Thus, 
the constraints given by 8 R = 0 are not independent 
of other constraints applied to 8 and the stepwise pro- 
cedure used to identify nonzero regional differentials 
need be applied to only R - 1 regions. As a result of 
the nonsingularity of G, it is readily apparent that if 
all differentials are zero, so that 8 = 0, then A[3 = 0 
and all like-process coefficients across regions are 
equal, implying the hybrid model is equivalent to the 
national model. 

Let 8 C refer to the collection of regional differentials 
determined by the stepwise procedure to be insignifi- 
cantly different from zero and let G c represent the 
rows of G corresponding to the K c = Y^f=i ^sf elements 
of d c , so that S c = G c p. The failure to reject the 
hypothesis that 8 C = 0 K c , where 0 K c is a .K^-element 
vector of zeros, leads to the relation 



G c p = 0 K c, 



(All) 



a condition that implies K 0 linear constraints across 
all the regional fixed-effects process coefficients. It is 
these constraints that are applied to the fixed-effects 
model to obtain the hybrid regionalized national 
model. 

A consistent estimate of the covariance matrix for 
the constrained model takes the form (Gallant, 1987; 
see the Supporting Information for a derivation) 



-l 



y(p) = V 0 - V 0 G c '(g c V 0 G c ') A G c V 0 , 



(A12) 



where Vo = a 1 
being the partia 



X>ln^ /; (/ 



derivative of the model estimate of 
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load for monitored reach i (of a total of N monitored 
reaches that supply the dependent variable load data 
used to estimate the SPAEROW model; see the 
description of the SPARROW model in the Support- 
ing Information) with respect to the parameter vec- 
tor, with derivatives evaluated at the constrained, 
WNLS estimates of the fixed-effects process coeffi- 
cients, p, and <t 2 = N^ 1 YaLi °f ■ VoP is positive defi- 
nite (positive semidefinite if any of the coefficient 
estimates are constrained according to model con- 
straints, such as that the source coefficients must be 
nonnegative), and Equation (A12) conforms with 
the linear constraints in Equation (All) in that 
premultiplying V(/j) by G c results in a zero matrix. 



SUPPORTING INFORMATION 



Additional Supporting Information may be found 
in the online version of this article: 

Data SI. Description of the SPARROW methodol- 
ogy, a derivation of the covariance matrix for the case 
where some of the coefficients are set to predeter- 
mined values, a derivation of the covariance matrix 
for the case where coefficients are estimated under 
linear constraints, and tables containing the esti- 
mates and related statistics of the TN and TP 
national, regional fixed-effects, and hybrid models, 
are available as part of the online article. 

Please note: Neither AWRA nor Wiley-Blackwell 
is responsible for the content or functionality of 
any supporting materials supplied by the authors. 
Any queries (other than missing material) should be 
directed to the corresponding author for the article. 
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